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ABSTRACT 


Solutions  to  a  nonlinear  wave  equation  were  analyzed  for 
their  stability.   The  wave  equation  is  a  Klein-Gordon  equa- 
tion with  the  mass  replaced  by  the  square  of  the  wave  func- 
tion.  This  wave  equation  has  propagating  solutions  which 
are  unbounded  or  periodic,  depending  on  the  sign  of  the  non- 
linear term  and  the  propagation  speed  which  can  be  sub-  or 
super-light  velocity.   The  stability  of  the  Deriodic  sub- 
light  velocity  solution  was  investigated  by  the  method  of 
characteristic  exponents  and  was  found  to  be  indifferent. 
Liapounoff's  direct  method  and  Sturrock's  analysis  of  the 
dispersion  relation  combined  with  a  WKB  technique  were  ap- 
plied to  a  linearized  perturbation  on  a  static  solution  of 

2 

the  field  equation.   The  periodic  solution  with  (3   <  1  is 

stable,  while  the  method  of  characteristic  exponents  gives 
indifference.   The  super-light  velocity  solutions  are  un- 
stable.  Due  to  the  limitations  of  the  approximations,  it 
could  not  be  determined  whether  the  instabilitv  is  absolute 
or  convective. 
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INTRODUCTION 

In  recent  years  there  has  been  increased  interest  in 
nonlinear  wave  equations  in  connection  with  theories  of  ele- 
mentary particles.   Such  equations  arise  when  one  attempts 
to  describe  a  particle  or  a  system  of  particles,  using  the 
theory  that  the  particle  itself  arises  from  self-interactions 
within  the  wave  function.   Heisenberg's  nonlinear  spinor 
theory,  [1]/  a  quantized  version  of  such  a  theory,  was  able 
to  reproduce  certain  features  of  the  spectrum  of  elementary 
particles . 

Nonlinear  differential  equations  have  solutions  that 
are  fundamentally  different  from  the  solutions  of  linear 
equations . [2]   They  can  disappear,  become  singular,  or  arise 
from  nothing.   A  central  problem  lies  in  the  stability  of 
the  solutions. 

The  purpose  of  this  investigation  is  to  study  a  possible 
method  of  stability  analysis  as  applied  to  the  solution  of  a 
wave  equation  which  may  serve  as  representative  of  a  whole 
class  of  nonlinear  wave  problems. 

Three  methods  of  stability  analysis  will  be  emoloved. 
Two  of  these,  Liapounoff's  method  and  the  method  of  charac- 
teristic exponents,  [3],  are  interrelated.   Liapounoff's 
method  of  stability  analysis  concerns  the  stability  of  the 
equilibrium  points.   A  theorem,  also  due  to  Liapounoff, 
states  that  the  results  regarding  stability  of  the  linear- 
ized theory  also  hold  for  the  nonlinear  case,  provided 
Liapounoff's  method  is  applicable  (the  characteristic 


equation  does  not  have  pure  imaginary  roots) .   The  method  of 
characteristic  exponents  concerns  the  stability  of  periodic 
solutions . 

The  third  method  is  the  analysis  of  the  dispersion  re- 
lation as  described  by  P.  A.  Sturrock  [4],  which  was  devel- 
oped from  stability  problems  in  plasma  physics.   This  method 
is  also  only  applicable  in  the  linearized  theory. 

The  representative  differential  equation  to  be  used  is 

V2$  _  ex2<J>5  -  ia  Ml  -  o  , 

c   9t2 

which  will  be  derived  in  Chapter  I  from  the  Klein-Gordon 
equation. 

The  particular  questions  being  asked  are: 

1.  What  are  the  static  solutions  of  the  equation 
for  given  symmetries? 

2.  Are  time  dependent  solutions  to  perturbations 
stable,  or  for  what  wave  lengths  are  they 
unstable? 

Chapter  I  of  the  paper  contains  the  derivation  of  the 
wave  equation  to  be  used,  an  analysis  of  the  types  of  solu- 
tions and  some  of  the  pertinent  solutions.   In  the  second 
chapter  selected  solutions  and  perturbations  are  analyzed 
for  stability.   Chapter  III  is  a  summary  of  results  and  con- 
clusions drawn  from  those  results.   Two  appendixes  attached 
to  the  paper  contain  expansions  or  alternatives  to  solutions 
derived  in  the  paper  itself. 


CHAPTER  I 
THE  NONLINEAR  WAVE  EQUATION 

1.1   DERIVATION  OF  THE  NONLINEAR  EQUATION 

The  wave  equation  chosen  as  a  starting  point  for  the 
derivation  is  the  Klein-Gordon  wave  equation  for  a  free 
particle : 

♦  '-  0  '.         (1) 


v*  -  ^t!  -  L-  ii 


c2  dt2 


I 


Two  reasons  may  be  advanced  for  this  choice.   First, 
the  Klein-Gordon  equation  is  the  basic  equation  for  relativ- 
istic  quantum  mechanics  reflecting  relativistic  energv  con- 
servation.  Therefore  a  nonlinear  equation  based  on  this 
ought  to  contain  similar  features  which  could  be  interoreted 
on  the  basis  of  the  known.   Second,  both  the  electromagnetic 
and  Yukawa  static  potentials  can  be  derived  from  the  relativ- 
istic equation  by  setting  the  mass  term  equal  to  zero  or  the 
mass  of  a  pi-meson,  respectively.   Thus,  one  can  hope  to  find 
in  the  static  solution  of  the  nonlinear  theory  the  analogues 
to  the  known  static  potentials. 

To  derive  a  nonlinear  wave  equation,  an  anproach  is 
taken  similar  to  the  one  that  led  Heisenberg  to  his  nonlinear 
spinor  equation.  [1] .   Instead  of  making  the  mass  term  a  con- 
stant, we  choose  to  write  m  ~  <$>2  ,  and  replace  the  entire  term 
by 

£\2<$>*         a 

where  e    =  ±1. 

This  choice  appears  to  be  plausible  by  recalling  that 
in  quantum  mechanics,  the  probability  of  finding  a  particle 


at  the  position  r  is  P(r)dr  =  <j>*<j>dr.   The  particle  can  be 
thought  of  as  being  smeared  out  over  soace  and  the  require- 
ment that 

I     <J>*<f>  dV  =  1 
J  v 

means  that  the  particle  remains  intact.   In  a  certain  sense 
then,  the  probability  of  finding  a  particle  at  the  position 
r  could  be  construed  to  indicate  its  mass  density  at  the  posi- 
tion r.   The  wave  equation  to  be  used  is  therefore 

V24>  -  £A2d>5  -  ia  ill  =  0  .  (2) 

c  at2 

£  =  1  corresponds  to  a  modified  Klein-Gordon  ecruation  with 
the  mass  replaced  by  self-interaction  of  the  field.   e  =  -1 
represents  the  equation  of  an  electromagnetic  potential  with 
a  source  term  given  by  the  field  self-interaction. 

1.2   BEHAVIOR  OF  THE  STATIC  SOLUTIONS 

The  static  potential  function  is  the  solution  of  the 
equation 

V2cf>0    =    eA2cf>*     .  (3) 

The  equation  can  be  solved  using  spherical  coordinates 
and  assuming  that  the  solution  is  spherically  symmetrical. 
The  equation  to  be  solved  is  therefore: 


2  d\ 

r   dr 


■•♦O-.xN? 


dr' 
To   proceed,    assume   the    solution   is   of   the    form 

<f>0    =    rnF(u) 


where  u  is  an  unknown  function  of  r.   By  substitution  one  ob- 
tains : 


»  SLl£  („.,*  +  r»  <»  u-  +  2(n+l)rn-1  £  u-  + 


r 


du2  du  ^ 

+  n(n+l)rn"2  F  =  eX2r5nF5   . 

One  may  choose  n  and  u  arbitrarily,  and  then  determine  F. 
If  u  =  In  r,  then  the  differential  equation  becomes 


n-2  d^F     n-2  dF    0/  ...  n-2  dF     ,  , , ,  n-2       ,  2  sn  5 
r r    —  +  2(n+l)r    —  +  n(n+l)r    F  =  eX2r5  F5. 

du' 


,2        du  du 


Choose  n  so  that  n  -  2  =  5n;  n  =  -%.   The  differential 
equation  then  reduces  to: 


%F  -  eX2F5  =  0  . 


d  F  ,2^5 


du" 


Let 


dF 

z  =  an 


dz    1  d2F 
dF  "  z 


or 


du' 


dz    d2F 

Z  dF  :  71 

du 


Therefore : 

z  ||  =  eX2F5  +  %F. 

z2     eX2F6    F2     . 

Note  that  this  equation  resembles  the  energy  integral  of  clas- 
sical mechanics  where  z  is  analogous  to  velocity,  the  function 
of  F  to  a  negative  potential  function,  in  which  F  plays  the 
role  of  the  space  coordinate,  and  h  to  the  total  energy.   One 


can  write  the  equation  in  the  form: 


?h+^AlFl  +%p2  u  (5) 

One  can  ask  what  a  plot  of  z  versus  F  (a  coordinate- 
velocity  phase  diagram  in  the  mechanical  analoq)  would  look 
like  for  various  values  of  h  and  e,    provided  one  requires  that 
z  and  F  remain  real.   A  qualitative  analysis  yields  the  re- 
sults in  the  following  sections. 

1.2.1   THE  CASE  e  =  1 
If  e  =  1,  then 


2p6 

z  =  ±  \/2h+"  £   +%F 


If  h  =  0,  then  for  large  F 

z  ~  ±  |F3 |  . 

If  h  >  0,  z  has  a  minimum  of  /2h  when  F  =  0,  but  other- 
wise behaves  much  in  the  same  manner  as  for  h  =  0. 

If  h  <  0,  z  can  take  on  all  values,  and  for  large  F  var- 
ies as  ±  |F3|  .   In  this  case,  however,  when  z  is  zero  the 
absolute  value  of  F  is  greater  than  zero. 

That  is 

| F  min |  =  a  >  0  . 

A  little  algebra  will  yield  the  one  real  value  of  a  as 

a  =  (A+B) s    , 

1 

A  =   |r-    [12b*A    +     Vl44b^X2+l]    3 

1 
B.=  2X    [12boX    "     Vl44bJX2+l]     3 

and  bn   =    I 2^ I •      T^e    results    are    shown   graphically 

in   Figure    1. 
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FIGURE  1 


,dF 


Phase  Diagram  (-=—  vs.  F)  for  Static  Case,  e  =  1 

On  each  curve  the  curve  parameter  is  u.   The  arrows  indicate 
the  direction  in  which  the  parameter,  u,  increases.   The  ac- 
tual values  of  the  parameter  along  the  curves  depend  on  the 
particular  solution. 

1.2.2   THE  CASE  e  =  -1 

In  this  case  equation  (5)  becomes 


2h-Al|l  +  %F2  . 


If  h  >  0,  a  phase  plot  will  yield  a  closed  looo  about  the 
origin,  where  values  of  z  and  F  outside  the  loop  are  forbidden, 
The  characteristics  of  the  loop  are  dependent  on  the  size  of 
h  and  A . 
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First,  one  notices  that 

I zl  =  /2h 


when 


F  =  0  , 


4X 


k 


If  z  is  set  to  zero,  one  may  solve  an  equivalent  cubic 

for  F  max, 

3  _  3v  _  0  _  0 


4A 


A' 


F  max  =  ±  (A+B)^  if  144A2bJ  >_  1  . 


A  =  ±J    l12X    bo  +  Vl44X2bJ  -  1]  3, 


B  =  iy  [12Ab2  -  Vl44A2b£  -  1]  T 


and 


b2  =    2h  . 


FIGURE  2 

dF 
Phase  Diagram  (-5—  vs.  F)  for  Static  Case,  e  =  -1,  h  _>  0 


When  144A2b£)  <  1, 


F  max  =  ±  [j   cos  9/3] 2    ,  cos  0  =  12b2A  . 
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There  is  only  one  root  for  F  since  if  144A2b£.  >  1,  two 
roots  of  v  are  imaginary;  if  144A2b!t    1,  two  roots  of  v  are 
negative  and  would  yield  an  imaginary  F.   These  results  are 
illustrated  in  Figure  2. 

If  h  =  0,  the  double  loop  in  Figure  2  is  obtained.   This 
case  is  only  a  degenerate  case  of  h  >  0.   The  two  nodes  at 

F  =  0  meet  at  the  origin.   In  this  case 

/  3   |  h 
F  max  =  ± 


IX2 
When  h  <  0,  equation  (5)  may  be  written 


2T?6 

z  =  ±"\/-b^  -'±-E—   +%F 


where 


bo  -  I 2h 


F  now  has  both  a  maximum  and  a  minimum.   The  absolute  value  of 
the  minima  must  be  greater  than  zero.   The  effect  is  to  sep- 
arate the  two  lobes  for  h  =  0  in  Figure  2  and  move  them  along 
the  F  axis.   From  the  solution  for  the  roots  of  the  cubic, 

3    3     ^  3b0     n 

it  can  be  seen  that  there  is  no  real  F  if 

144A2b£  >  1 

because  the  roots  for  v  are  either  negative  or  imaginary.   If 

144X2bJ  <  1  , 
there  will  be  two  solutions  for  F.   These  represent  the  maxima 
and  minima  for  F.  -. 

F  max,  min  =  ±  (jcos    6/3)*  ,  ±  Ijji^l   sin  9/3-cos  9/3) ]3 

where  cos  6  =  -12Ab2  . 
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If  144A2b!l  =  1,  there  is  one  real  root  for  F, 

This  represents  a  point  solution  for  F.   The  results  are  di- 
gramed  in  Figure  3. 


h<0 


k 


:  /i_ 

2\ 


FIGURE  3 


.dF 


Phase  Diagram  (-^—  vs.  F)  for  Static  Case  z   =  -1,  h  <  0 


All  the  solutions  for  <$>Q   will  be  periodic  except  when  h  =  0 , 
or  when  F  recedes  to  a  point  when  h  <  0. 


1.2.3   ANALYTIC  SOLUTIONS  FOR  <j>-  WHEN  h  =  0 

The  two  most  interesting  cases  in  the  last  two  sections 
are  those  for  which  h  =  0. 

When  e  =  1 ,  h  =  0 ,  one  obtains  the  equation 

(^2    -  hill   +  ^F2 
The  equation  may  be  solved  by  letting 


F  =  v 


-h 


Then 


3  0   -  3v*  +  4X 
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dv        .     , 
7Y2-7:  =  In  r/< 


.  -1  v/3    ,    , 
sinh    — 7TT-  =  In  r/x  , 


2X 


-1 


which  upon  transposing  the  sinh    and  using  the  exponential 
transformation  for  the  sinh,  one  obtains : 

X(r2-<2) 


v  = 


*0  = 


/3<r 


/Ik 


A(r2-<2). 


h 


This  solution,  (see  Figure  4),  represents  a  potential  of 
infinite  height  at  r  =  < .   Suppose  this  potential  represented 
the  potential  from  a  particle.   Such  a  potential  would  indicate 
that  the  particle  consisted  of  a  central  core  of  radius  k  in 
which  there  is  no  real  potential.   Outside  the  core,  the  po- 
tential varies  much  as  a  classical  1/r  potential.   The  concept 
of  a  hard  core  is  not  new.   Arguments  have  been  advanced  that 
such  a  core  might  exist  in  nuclei  between  particles.  [5] 


* 


0 


< 


FIGURE  4 


Static  Solution  (4>Q  vs.  r)  e 


=1,  h=0,  r>K 
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It  is  surprising  that  one  gets  a  long  range  potential  in- 
stead of  the  short  range  Yukawa  potential  of  the  Klein-Gordon 
equation. 

Using  the  solution  for  ((>-.  one  can  determine  the  applicable 
regions  of  the  phase  diagram  in  Figure  1. 

Since 


F  =  V  2 


F  =  ± 


/3k: 


X(r2-K2) 

When  F  =  0,  r  =  °°.   That  is,  the  saddle  point  in  Figure  1 

belongs  to  u  =  °°. 

When  considering  the  derivatives,  one  obtains 

dF        1  dF 
du       r  dr 

dF 
When  F  >  0,  -^—  <  0.   Therefore,  the  applicable  region  o* 

dF 
the  phase  diagram  is  the  fourth  quadrant.   When  F  <  0,  -5—  >  0. 

The  applicable  region  in  the  phase  diagram  is  the  second  quad- 
rant. 

The  first  and  third  quadrants  must  then  belong  to  a  dif- 
ferent solution.   In  arriving  at  the  solution,  the  equation, 

3(5u")2  =  3v'  +  4X'  ' 

was  derived  and  the  solution  was  obtained  using  the  positive 

root  for  the  derivative.   If  now  the  negative  root  is  used, 

one  will  eventually  get 

.,  -1  v/3    ,    , 
sihn   — 2Y  =  In  K/r  . 

The  solution  for  <}>n  in  this  case  is 


*0  =  ± 


/3k 


A(K2-r2) 
16 


h 


This  solution  represents  a  potential  well,  the  general  features 
of  which  are  shown  in  Figure  5. 


* 


0 


FIGURE  5 


Static   Solution    (<f>Q    vs.    r)    e    =    1,    h 


=  0 ,  r  <  k 


This  solution  is  represented  in  Figure  1  by  the  first  and 
third  quadrants.   The  saddle  point  belongs  to  r  =  0 ,  or  u  =  -°° 

If  one  solves  the  static  equation  for  e  =  -1,  h  =  0, 
quite  a  different  result  occurs. 


Let  F  =  v 


Then 


3(f*V)    =   3v2   _  4A2 

du 


which  has  the  solution 

,  -1  v/3    ,    . 
cosh   ~~TT  ~   ln  r/ic  » 

where  k  is  the  integration  constant.   Transposing  the  inverse 
cosh,  substituting  the  exponential  equivalent  for  the  cosh,  and 
then  solving  for  v,  will  yield 

A(r2+K2) 


v  = 


x/3 


<j>0   =    (rv) 


-h 


/3 


Ur2  +  <2) 


l2 


7) 
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* 

Equation  (7)  represents   an  equation  which  bears  a  strik- 
ing resemblance  to  a  static  electromagnetic  field  potential. 
Besides  the  resemblance  to  coulomb's  law,  this  solution  has 
another  point  of  similarity  with  electric  potential  fields. 
The  differential  equation  from  which  it  was  developed  is  equa- 
tion (3)  with  e  =  -1, 

V2(f)0  =  -A2^  (8) 

One  might  compare  this  with  the  static  potential  equa- 
tion with  a  source  term, 

V2*  1    -  £-   , 

Eo 

where  cj>  represents  potential,  p  the  charge  density,  and  en  the 
permittivity  of  free  space.   This  potential  function,  and  the 
wave  equation  it  implies,  shall  be  used  for  the  time  dependent 
stability  analysis  in  the  following  sections.   Other  solutions 
of  the  static  wave  equation  may  be  found  in  Appendix  A. 

Assume  that  the  <j>fl  of  equation  (7)  represents  the  coulomb 
potential  from  an  electron.  One  might  ask  what  would  the  con- 
stants k  and  X   be,  and  what  would  they  represent. 

To  determine  the  constants  two  relations  will  be  required, 

j 
The  first  is  derived  from  the  fact  that  in  electromagnetism 

P 


and 


where  q  is  charge.   In  this  case  q  =  e,  the  electronic  charge. 

*  Chandrasekhar  [6]  quotes  this  solution  in  the  theorv 
of  stellar  structure. 
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But  for  the  nonlinear  wave  equation, 


v2*o  =  -*2*o  =  h 


o 

and   in    spherical    coordinates 

r7T  r2l\ 


</3 


r2+<2 


—  =      de 


o 


J  o 


r 


5      5 

1         2    ^T  r2    dr   sin    9 
~T   K       3      k 

(r  +k  ) 


dip 


0 


5      5 

4tt  •  3      k 
A* 


Let  r  =  k  tan0. 
Then  the  integral  becomes 

-tt/2 


sin 


r  °° 


r   dr 


(r2+K2) 


cos  e  de 


e_ 

E, 


and 


^  A- 
4ttk*3* 

A^ 


% 


=  9'1094fTe 


(£)   =  34-109e   . 

Conversations  with  those  performing  the  collidinq  beam  ex- 
periment at  Stanford  University  indicate  that  coulomb's  law 

-17 
holds  to  within  2%  at  r  =  3*10    meters.   One  may  use  this 

data  to  obtain  a  second  relation  between  k  and  A. 

From  the  ratio  for  k/A   evaluated  above, 

.     9«109e 

*°  =  (r2+K2)^    ' 
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Now,  let  £,  and  £?  be  defined  so  that 


?1    r  '  ?2 


T 


U2  +  K2)S 

A  plot  of  these  two  functions  is  shown  below  in  Figure  6 


r0=3.10 


-17 


FIGURE    6 
Comparison   of   Coulomb   Potential    and   <f>0    For   e   =   -1,    h  =    0 


At  rfi ,  Co  cannot  be  more  than  2%  below  £,  •   That  is, 
?1  >  S2  >  .98^  at  rQ  . 

This,  in  effect,  puts  an  upper  bound  on  k. 

Let 

C2  =  .98C]_  at  rQ. 


Then: 


|i  •  10" 


[9-10  "  +  k2]  s 


=  3T~ 
-17 


,-17 


k  max  =  .6  10 
The  range  of  k  is 

0    k    .6*10 

using  the  upper  limit  of  k  to  obtain  X,  one  would  get 

X   =  5. 
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1.3   BEHAVIOR  OF  THE  SOLUTIONS  OF  THE  TIME  DEPENDENT  EQUATION 
One  can  find  propagating  solutions  to  the  time  dependent 
equation  in  one  dimension.   Such  a  solution  will  give  some 
idea  of  the  behavior  of  the  solutions  for  more  complicated 
symmetries.   The  time  dependent  wave  equation  in  one  dimen- 
sion is: 

9x2    c2  9t2 
Let 

T  =  Ct  , 

S  =  X  -  3-T   , 

»-5 

Assume    that    <f>    is   of   the    form 

0    =    <J>(S)     , 
i.e.,    we    confine   our   attention   to   propagating    solutions. 

Then 

92(fr   _   d2(j) 

and 


9x2 

dS2 

92cj) 
9x2 

g2d2* 
dS2 

The 

dif ferenti 

al 

equation    is 

d2(J) 

eA2c}>5 

Let 


Then 


dS2  1-6 


d<f>         , 
dt  =   *    • 


$±       I  111 
d<D        *  ds2 
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(9) 


and 


*  d*      1-B2 


As  in  section  1.2  the  first  integral  yields  an  equation 
similar  to  the  energy  equation  in  classical  mechanics, 


6(l-32) 


=  h    , 


or 


iff   =    ±  .  /2h+ 


eA2<}>6 


V  3(l-32) 

A  qualitative  analysis  of  this  equation  using  the  methods 
in  section  1.2  is  presented  in  the  following  section. 


1.3.1   DISCUSSION  OF  \\)    vs.  <J> 

If  in  the  equation  for  ip , 

e  =  1,  and  3 2  <  1  , 
curves  of  cf)  vs.  ip  will  be  the  same  as  those  in  Figure  1  with 
F  replaced  by  <j>  and  z  by  iJj.   Again  this  case  represents  solu- 
tions which  are  unbounded. 

If  e  =  1,  32  >1,  the  second  term  in  the  equation  for  ip 
is  negative  and  h  must  be  greater  than  0.   The  phase  diagram 
will  look  something  like  Figure  7. 


FIGURE    7 


.dcj> 


Phase   Diagram    (-—  vs.    <j))     for  Wave   Equation  With 
e  =   1,    $2   >   1,   or  e  =   -1,    32    <   1 
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One  has  the  result  then,  that  for  e    =    1,  real  sub-liqht 
velocity  solutions  exist  which  are  unbounded.   The  only  peri- 
odic solutions  are  those  which  exist  for  suoer-light  veloc- 
ities . 

No  other  cases  need  be  considered,  since  if  z   =  -1,  one 
merely  needs  to  invert  32  in  the  last  two  cases  to  obtain  the 
same  results.   That  is 

62  >  1,  £  =  -1   corresponds  to   62  <  1,  £  =  1  , 
and     B2  <  1,  e  =  -1   corresponds  to   62  >  1,  e  =  1  . 

1.3.2   ANALYTIC  SOLUTIONS  OF  THE  ONE  DIMENSIONAL  EQUATION 

Because  of  the  redundancy  of  cases  the  solutions  shall  be 

carried  out  with  32  <  1  . 
Let 

2h   =    6b2     ;    6    =    ±1    . 


d£ 
dS 


=   bn-\/6+   £i^ =   bn[6  +  £a^6]%    ,  (10) 

0V  3(l-62)b2  °  ° 


where 


a~    = 


Let 


then 


Let 


'3(1-B2)b2 


£   =    a04>    , 


1  =    a0bQVT^F 


then 


dv 


=%a0b0(S-S0) 

[v(6+evd) ] 
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where  the  limits  on  the  integral  have  been  left  open  since 
they  depend  on  e . 

If  e  =  <5  =  1, 
the  solution  is 

cos  9  =  cn(Yfk)  , 
where  en  is  an  elliptic  function  with  modulus  k.   The  variables 

are 

1  +  (l-/3)v 

cos  9  =  - ■ —    , 

1  +  (l+/3)v 


2  +  /3 


kz  = 


% 


Y  =  2a0bQ34(S-S0)  . 


Substituting  for  cos  9  in  the  solution,  one  obtains  v  as 

1-cn  y 


v  = 


and 


1   %    1 

4)  =  —  v  =  — 

a0      a0 


(/3-l)  +  (/3+l)cn  y 
1-cn  y 


(/3-l)+(/3+l)cn  y. 


k 


Let   6  =  -1,  e  =  1. 


Then 


-v 


dv 


Substitute   v  =  —  . 


1  [v(v3-l)]% 


1 


2a0b0(S-S0) 


w 


Then 


■w 


dw 


(l-w4)% 


2a0bQ(S-S0) 
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The  solution  is 

w  =  en [2/2aQb0 (S-SQ) ]  =  en  y 

with  modulus 

k2  =  % 

1 


<f>  = 


a0/cn~Y 


One  other  solution  may  be  obtained  when  b2   or  h  =  0. 

Then  must  be  greater  than  zero.   In  this  case  the  equa- 

1-32 

tion  becomes : 

dd)       Ad)  ,  3 

dS    [3|1-32M 


■? 


where  a  can  be  either  positive  or  negative. 

<j>  =  ±[2a(SQ-S)]"%  . 

These  solutions  are  all  unbounded.   It  might  be  of  inter- 
est to  construct  a  phase  diagram  similar  to  that  of  Figure  1 
with  the  limiting  values  of  S  indicated.   The  h  =  0  case  is 
represented  by  the  last  solution  obtained.   When  a  is  positive, 
and  cf>  is  negative,  the  third  quadrant  curve  is  obtained.   The 
saddle  point  belongs  to  S  =  -00  in  this  case.   If  a   is  nega- 
tive, the  second  and  fourth  quadrant  curves  are  applicable . 
The  saddle  point  belongs  to  S  =  °°.   The  results  are  shown  in 
Figure  8.   The  arrows  indicate  the  direction  of  increasing  S. 

When  h  ^    0,  the  other  curves  are  obtained,  in  general, 
using  the  same  methods  as  in  section  1.2.1.   Ifh  >  0,  •—   is 
a  minimum  when  <j>  =  0 .   This  occurs  when  en  y  =  1  or  S  =  SQ . 
The  derivative  approaches  infinity  when  <J>  approaches  infinity 
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FIGURE  8 
Phase  Diagram  (-—   vs.  <J>)  for  Wave  Equation  with  e=  +1,  32  <  1 


or  when 


/3-1  =  (/3+l)cn  y  / 


for  the  solution  obtained.   The  en  y    lies  somewhere  in  its 
second  quarter  cycle.   That  is 

K  <  y  <  2K  , 

where  K  is  the  complete  elliptic  integral  of  the  first  kind. 

S  would  have  some  corresponding  value,  S, ,  as  labeled  on  ^ig- 

def) 


ure  8 .   When  h  <  0 , 


dS 


0  when  E,   =  1.   This  corresnonds  to 


en  y  =  1/  or  S  =  S..   The  derivative  approaches  infinity  when 

en  y  =  0  . 

This  occurs  when  en  y  is  at  its  quarter  evele,  or  S  =  Sj 
This  is  the  label  used  in  Figure  8. 

These  solutions  became  infinite  for  certain  values  of  S. 
Adjacent  to  these  infinities  the  function  is  not  real.  It  is 
therefore  not  possible  to  construct  at  any  time,  t,    a 
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perturbation  which  is  real  for  all  space.   These  solutions  do 
not  allow  a  ready  physical  interpretation. 
If6=l,  e=-l,  the  equation  becomes 


dv 


[v(l-vJ)] 


— ¥=  2a0b0(S-S0)  . 


The  solution  is 


■% 


cos  0  =  cn(2aQbQ3  (S-S.)  =  en  y    , 
where  the  modulus  of  the  elliptic  function  is 

k2  = 
and 


2-/3 


cos 


1 —  ' 

1- (14-/3)  y 

l+(/3-l)v 


Upon  making  the  proper  substitutions  and  solving  for  v, 
one  obtains  <j> ; 

(11) 


,         1    [  1-cn   y 

ao|l+/3+(/3-l)cn 


Y 


This  solution  represents  the  only  type  of  solution  that 
is  oscillatory  in  nature. 
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CHAPTER  II 
STABILITY  OF  SOLUTIONS  TO  THE  WAVE  EQUATION 

In  this  section,  the  stability  of  the  solutions  to  the 
wave  equation  is  investigated.   The  time  dependent  periodic 
solution  is  investigated  using  the  method  of  characteristic 
exponents.   The  following  sections  will  discuss  the  stability 
of  wave-like  perturbations  on  the  static  solution  for  e  =  -1. 
The  methods  to  be  used  are  the  characteristic  exponent  method, 
Liapounoff's  method  and  Sturrock's  stability  analysis.   Each 
of  these  methods  shall  be  described  as  they  occur. 

2.1  THE  UNBOUNDED  TIME  DEPENDENT  SOLUTIONS 

The  solutions  obtained  in  section  1.3.2  with  e.   =  1  in 
equation  10  all  yield  unbounded  functions,  not  real  for  all 
r  and  t.   These  solutions,  given  in  section  1.3.2,  are  all 
basically  unstable  and  shall  not  be  discussed  further. 

2.2  STABILITY  OF  THE  PERIODIC  SOLUTION 

The  solution  of  equation  10  in  section  1.3.2  yielded  a 
periodic  solution  (equation  11)  when  e  =  -1.   The  solution 
shall  be  investigated  for  stability  using  the  method  of  charac- 
teristic exponents. 

Suppose  that  one  has  a  system  such  that 

^  =  P(u,v)  , 


and 

dv 


dt  =  Q<u,v)  . 
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Suppose  the  solutions  to  these  equations  are  known  to  be 

ul  =  ui^t^  and  vi  =  v,  (t). 

Since  these  solutions  are  assumed  to  be  periodic,  they 
must  have  a  period  which  shall  be  called  T. 
Then  define 

h   -  1   f  r8p    3QiH4- 

h  "  f  [TZ+   9v"]dt  • 


0 

Then  the  periodic  motion  is  stable  if  h  >  0,  and  unstable 
if  h  <  0  [4].   Applying  this  approach  to  the  specific  problem 
at  hand  one  obtains 

i  =  *  =  p  - 

and 

2  a  5 


St   =  -  X    *       =    Q 

1-6" 


dS      ,  n2 


These  differential  equations  are  essentially  equivalent  to 
derivatives  with  respect  to  time  if  the  analysis  is  carried 
out  for  constant  position,   that  is, 

it  =     I  3i 

dS    "  3  dx 
at  a  given  value  of  x  since 

dS  =  dx  -  3  dx  . 
Carrying  through  the  other  operations  one  finds  that 


and 


1£  -  11  =  n 


3Q  =  a  (Zilil)  .  o  . 


8  ip    9i|; 


i-e2 


Therefore  h    0. 
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The  theory  then  predicts  that  the  motion  is  indifferent 
to  perturbations.   That  is,  the  perturbation  results  in  a  new 
amplitude  of  oscillation  which  will  be  maintained. 

2.3   STABILITY  OF  THE  PERTURBED  STATIC  FUNCTION 
2.3.1   THE  WAVE  EQUATION  FOR  THE  PERTURBATION 

One  can  derive  a  wave  equation  for  a  perturbation  applied 
to  the  static  solution.   As  stated  in  section  1.2.3,  only  the 
case  for  e   =   -1  shall  be  analyzed  so  the  static  ootential  is 
represented  by  equation  (7). 

The  wave  equation  of  interest  is 

V24>  +  A2cJ>5  -  ia  i!±  -  0  .  (12) 

C   3t2 

Assume  that 

where  ip  is  very  small. 
Since 

equation  (12)  becomes 

v24)n  +  v2^  +  A2(<j)n+i|,)5  -  i2  i!l  =  o. 

0  0        c   9t2 

The  third  term  may  be  expanded  to 

X2  (cj)S  +  5<i>^+.  .  .+ijj5) 

But   if   \p    is    very    small   one  may   neglect   all   but    the    first   two 
terms.      Then,    since 

V2^    +    A2c^    =    0    , 


<J>    =    <t>0    +    </>    / 


3t  ' 


one   may   write 

V2^   +    5A2(j)>    -   i2 

at 


V2^   +    5A2(j)>    -   ^   111  -    0.  (13) 

U  C  r\    j_    2 
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Equation  (13)  shall  be  used  as  the  wave  equation  for  a 
perturbation  applied  to  <J>n  . 

This  problem  is  of  interest  because,  if  <t>0  represents  any- 
thing like  a  particle,  then  the  question  arises  whether  the 
nonlinear  self-interaction  allows  a  stable  particle  or  whether 
it  tends  to  break  up  because  of  application  of  perturbations. 

2.3.2   LIAPOUNOFF'S  STABILITY  ANALYSIS  AND  STABILITY  OF 
EQUILIBRIUM 
Suppose  that  one  picks  a  point  in  space  and  proceeds  to 
analyze  the  reaction  of  a  perturbation  approaching  that  point. 
Then  <|>n  in  equation  (13)  is  a  constant.   Using  the  substitu- 
tions and  the  method  outlined  for  the  periodic  solution  in 
section  2.2,  one  obtains  for  a  one  dimensional  perturbation 

ill  =  .  !Al!i!  ,  <14> 

dS2       1-32 
where 

5X2^ 
1-32 
is  a  constant. 

Since  the  equation  was  derived  with  x  as  a  constant,  the 
derivative  is  essentially  with  respect  to  time.   One  can  now 
do  a  local  stability  analysis.   Where  the  word  local  is  used 
to  emphasize  the  fact  that  only  one  point  in  space  is  being  in- 
vestigated.  The  equation  derived  above  obviously  represents 
periodic  rti     n  if :  32  <  1,  and  an  unbounded  solution  if  B 

Howevea   a  nonlinear  wave  equation  has  been  converted  to 
a  linear  equation  by  assuming  that  \p    is  small.   One  might  ask 
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if  such  a  prediction  has  any  validity  for  the  nonlinear  problem, 
Liapounoff's  stability  analysis  will  yield  the  answer  to  this 
question. 


Let 


and 


and 


du    _  ,         N 
gt  =  Plu,v)  , 


g=  Q(u.v)  , 


dv   P    cu+dv+Q' (u,v) 
du  ~  Q    au+bv+P1 (u,v)  ' 

where  Q'  and  P1  have  no  first  order  and  no  constant  terms.   If 

u  and  v  are  near  zero,  the  last  equation  can  be  reduced  to 

dv   cu+dv 
du    au+bv  ' 

and  the  equation  has  been  linearized. 

Define  s  so  that 

a-s   c 

=  0  .  (15) 

b    d-s 

Liapounoff's  stability  criteria  are  as  follows,  if  the 

roots  s,  and  s-  of  the  determinant: 

(1)  are  complex  or  real,  and  Re(s,,  s-)  are  negative, 
then  the  solution  is  stable. 

(2)  are  complex  or  real,  and  at  least  one  of  the 
roots  has  a  positive  real  part,  the  solution  is  unstable. 

(3)  are  imaginary  and  Re(s,,  s~)  =  0,  then  the 
stability  analysis  is  questionable.   The  question  as  to  the 
solution's  stability  is  not  conclusively  answered. 

Applying  this  analysis  to  equation  (14)  one  proceeds 
as  follows: 
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Let 


dip 
dS 


=    S    , 


then 


d^    -^2»^ 


dS 


1-3 


=     OL\l>      , 


and 


d£        aip 

Then  in  equation  (15) 

c  =  a,  b  =  1,  and  d  =  a  =  0. 
Therefore  s2-a=0,s=±/a. 
If  B2  >  1,  a  >  0,  and  the  roots  of  s  are  real, 


s  =  ± 


5x2*S 


1-6 


* 


Therefore,  one  can  say  that  if  82  >  1,  the  solution,  as 
already  stated,  is  unstable. 

If  32  <  1,  the  roots  of  s  are  imaginary.   Therefore,  the 
analysis  does  not  give  reliable  information  concerning  the 
stability  of  the  motion.   However,  the  motion  is  periodic  and 
the  method  of  characteristic  exponents  in  section  2.2  may  be 
used. 

Applying  this  method  once  again 


d^ 


dS 


=  P(if>,£)  =  K    , 


2  A  ■+ 


^=  Q  = 
dS    g 


1-32 


Then 


8P    90    n 

9ip  "  8£  " 
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Therefore  m 


3P    9Q 

3ip    9? 


dt  =  0  . 


0 

Since  h  is  zero,  one  concludes  that  the  motion  is  indifferent 
to  perturbations.   This  result  does  not  change  if  we  take  into 
account  all  higher  order  terms  in  the  expansion  of  (<f>0  +  ty)  5 , 
still  considering  <j>0  constant.   The  question  arises  at  this 
point  as  to  whether  this  approximation  to  "local"  stability 
or  instability,  by  keeping  <j>n  a  constant  and  not  really  con- 
sidering the  perturbation  as  propagating  into  regions  of 
changing  <j>_,  might  invalidate  the  conclusions.   It  certainly 
is  not  giving  information  about  propagating  instabilities, 
known  also  as  convective  instability.   In  order  to  investigate 
this  point  the  method  in  the  next  section  will  allow  removal 
of  the  restriction  that  <J>0  remain  constant.   The  stability 
analysis  technique  to  be  used  is  well  known  in  plasma  physics. 

2.4   STABILITY  ANALYSIS  WITH  THE  WKB  METHOD 

Another  general  approach  to  stability  has  been  developed 
by  Sturrock.   It  requires  that  the  dispersion  relation  to  the 
wave  equation  be  obtained.   For  this  analysis  a  three  dimen- 
sional, spherically  symmetric  solution  shall  be  obtained  for  a 
perturbation  on  the  static  potential  function. 

2.4.1   AN  APPROXIMATE  SOLUTION  TO  THE  LINEARIZED  WAVE  EQUATION 
In  section  2.3.1,  the  linearized  wave  equation  for  a  oer- 
turbation  on  the  static  potential  function  was  obtained: 

V2ip  +  5A20  -  i2  i!i  -  0  . 
0     C   3t2 
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If  the  static  potential  function  (equation  7)  is  put  in 
this  equation,  one  obtains 

1  d2\p 


V2*  +  ilf  * 


2  \  2 


=  0   . 


(r^+Kl 


at' 


It  will  be  assumed  that  the  perturbation,  \p ,  is  spher- 
ically symmetric,  thus  simplifying  the  mathematics  involved, 
and  amounts  to  considering  s-states  only. 

The  equation  then  becomes 

ill   +  1  ii+  15k2i|>     _  1   8^p  _  Q 
9r2    r  8r    (r2+K2)2    °2  8t2 

To  solve  this  equation  we  assume  the  solution  to  be  of  the 

form 

\p   =  T(t)R(r)  , 

and  a  separation  constant  a2.   Then  we  may  write: 

1_  dJ_R    2_  dR      15<2      __1   d^T  _    2 
.2  +  Rr  dr    ,  j.  ?v!    Tc2 


R   dr 


(r'  +  KM 


dt 


d2T 
dt2 


2^2 


+  Ta'cr  =  0 


m     ±ioat 

T  =  e     ,  oo  =  ac  . 


d2R    2  dR    „  f    15 


+  — 


+  R 


dr: 


Let 


r  dr     I  (r2+<2)2 


+  a 


=  0 


Let   rR  =  u  ,   then 


d^u 


+  u 


15 


+  a' 


=  0 


dr' 


I (rz+Kz) 


2  \  2 


2,„2 


r  =  kx,  a  k   =  U      , 
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then 


dzu 

dx2 


+  u 


15 
.  (x2+l) 


+  Q1 


=  0 


The  exact  solution  to  this  equation  is  rather  difficult 
to  obtain.   Therefore,  the  WKB  approximation  is  used  to  get 
a  solution  which  may  be  used  for  stability  analysis. 

If  one  has  a  differential  equation, 

+  K2u  =  0  , 


d2u 


dx' 


then  the  WKB  approximation  is  to  use 


u  =  K-*ei/Kdx 


For  this  case, 


K2  = 


+  ^2  , 


and 


K 


■k 


(x2+l)2 


(x2+l)% 


2/„2,n2i? 


[1-5+0'  (xz  +  l)  z] 


k 


The  evaluation  of  the  integral  in  the  exponent  will 
yield  elliptic  functions  of  all  three  types.   Such  a  func- 
tion will  be  very  difficult  to  analyze  for  stability.   There- 
fore, this  integral  will  also  be  approximated.   To  do  this 
the  Q,2  ,  x  plane  shall  be  divided  into  two  regions.   One  near 
the  origin  and  one  further  from  the  origin.   (See  Figure  9) . 
(a)   Assume  that 

— ^ <  1  , 

fi2(x2  +  l)2 

i.e.,  Q,2  (x2+l)  2   is  large. 
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Then 


K 


15 
I (x2+l) 


K 


+    ft' 


k 


=  ft 


n 


Kdx    -    fix   + 


1    +    % 


15x 


4ft(x2+l) 


-12 +  i 

Ift2 (x2+i) 2 

15 
ft2(x2+l)2 


15    .       -1 

+  m  tan 


% 


(b)       Now   assume    that 


ft2 (x2+l) 2 
15 


<  1  , 


Then 


i.e.,  ft  (x  +1)   is  very  small. 


K  = 


x2+l 


i  +  ft 


2   (X2+1) 


15 


% 


Kdx  a  /l5  tan  Xx  + 


ft2   ,1„ 


2/15 


(i-xd+x)  . 


To  summarize: 
if 


15 


,2  /„2 


<  1   , 


JT  (x'+l) 


u  = 


(x2+l)% 


[15+ft2 (x2+l) 2] « 


EXP  i 


fix  + 


15x     15  ^   -1 
+  —  tan   x 


4ft(x2+l)    4ft 


if 


ft2 (x2+l) 
15 


<  1  , 


u  = 


(x2+l) 


[15+ft2 (x2+l2]% 


EXP  i 


/15  tan   x  + 


2/15 


(ix3+x) 


Quite  a  few  approximations  have  been  made  in  obtaining 
these  solutions.   Therefore,  the  regions  in  which  the  aoorox- 
imations  are  useable  must  be  determined.   The  condition  for 
the  validity  of  the  WKB  approximation  shall  be  derived  first, 
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x  axis 

Regions  of  Validity  of  WKB  Approximations 
Figure  9 
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In  the  differential  equation, 


d2u 
dx2 


+  K2u  =  0  , 


make  the  independent  variable  substitution 

y  =  /Kdx   . 
Then 

HHk2  +  K'^  +  K2u  =  0  . 

dy2       dy 


where 


d^u   KJ_  du  +  u 

.2    „  2  dy 


=  0  , 


dy 


K' 


K'  = 


dK 
dx  * 


Now  let 


u  =  vs,  s  =  K  2ely 


ds    1  ds 
dy  :  K  dx 


ds 
dy 


1 
K 


-% 


3 

K2 


+  iK% 


ly 


d2s 
dy2 


1  |  5  (K1 )  2      K" 
K  I  T     7     *  ~5 


K' 


^  +  2i^X  +  v 
dy2      d* 


K 

3  (K')2 

4  K^ 


iK' 
3 
2 


-  K 


k 


K 


-  %• 


IT 

K3 


ly 


=  o  . 


If  the  last  term  is  zero,  then  the  general  solution  for 

u  is 

u  =  K-%  (ae"1*  +  be1*)  . 

Therefore  the  particular  solution  for  u  assumed  by  the  WKB 
approximation  is  good  when  the  absolute  value  of  the  last 
coefficient  is  very  small. 
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If  e  is  defined  to  be  equal  to  this  term,  then  by  dif- 
ferentiating and  by  algebraic  manipulation  one  obtains  ^or 
this  expression: 

15[15+ft2 (x2+l) 2]  -  75ft2x2 (x2+l) 2 
[15+fi2(x2+l)2]3 

The  shaded  area  near  the  bottom  of  Figure  9  is  the  re- 
gion of  space  for  which  |e|  >  0.1  .   Within  this  region  the 
WKB  approximation  is  either  not  valid  or  of  questionable 
accuracy. 

Next,  consider  the  approximation  of  /Kdx.   The  approx- 
imations  used  are  those  for  the  square  root  of  (z+1) 2    . 

(z+1)  2   =  1  +  %z  -  %z2  ...  ,  z  <  1 

where  only  the  first  two  terms  were  used  in  the  expansion. 

The  error  is  not  larger  than  the  absolute  value  of  the 
third  term.   Again,  the  value,  0.1,  is  chosen.   For  case  a, 

h  z2  = ISl —  .    .i      . 

4ft2  (x2+l) h   ~ 
Making    the   expression   an   equality,    we    can  write 
ft2(x2+l)2    =    23.75    ,    or   roughly    , 

ft2    =   -Jtt_         . 

(x2+l)2 

If  the  same  procedure  is  followed  for  case  b,  one  obtains 

ft2  =_^L-   < 

(x2+l)2 
These  two  functions  are  plotted  in  Fiqure  9.   The  shaded 
region  between  them  represents  the  region  where  the  approx- 
imation to  the  integral  is  unreliable.   The  region  near  the 
origin  is  represented  by  case  b;  the  region  above  the  shaded 
area  by  case  a.  4Q 


2.4„2   STURROCK'S  STABILITY  ANALYSIS 

Suppose  that  the  solution  of  a  wave  equation  yields 


,     i(kx-oot) 

<P  =  e  / 


and 


k  =  k(w)  , 
a)  =  co  (k)  . 
If  a)  is  plotted  in  the  co .  ,  to  ,  plane  for  all  real  k,  then 
the  shape  of  the  plot  will  yield  information  about  the  stabil- 
ity of  the  wave.   Furthermore,  a  plot  of  k  in  the  k. ,  k   Diane 
2  r  i '   r 

for  all  real  go  will  yield  information  concerning  the  tvae  of 
instability  of  the  wave. 

For  purposes  of  clarity  the  possible  alternatives  shall 
be  divided  into  four  cases. 

Case  (a)  -evanescent  waves 

These  waves  are  stable.   The  system  is  stable  if 
for  real  k,  id  is  always  real.   This  situation  is  diagrammed 

in  Figure  10. 

go  . 

l 


-€E 


^= 


-<- 


to 
r 


FIGURE  10 

Dispersion  Relation  in  the  Complex  go  Plane 
for  Evanescent  Wave 


41 


Case  (b)  -  non-convective  instability 

This  case  represents  waves  in  which  a  disturbance 
grows  and  expands  about  the  point  at  which  it  occurs ,  but  does 
not  travel. 

Such  a  case  occurs  in  the  type  of  plots  diagrammed  in 
Figures  11  and  12.   Note  that  w./io   -+  0  as  k  -►  ±  «>,  and  that 
k  is  real  for  all  real  oo. 


FIGURE  11 
Dispersion  Relation  in  the  Comnlex  oo  Plane 
for  a  Non-convective  Instability 
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FIGURE  12 
Dispersion  Relation  in  the  Complex  k  Plane 
for  a  Non-convective  Instability 
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Case  (c)  -  convective  instability 

This  case  represents  a  wave  disturbance  which  grows 

and  propagates.   The  to.  ,  to   plot  is  the  same  as  Figure  11, 

but  the  plot  now  also  has  complex  values  of  k  for  real  to  and 

k./k   -►  0  as  to  -*  ±  °°.   The  plot  of  k.  ,  k  would  be  the  same 
l'  r  1    r 

as  Figure  11  with  the  appropriate  changes  in  labeling  of 

axes . 

Case    (d)       In   this    case   to .  /to      -*■   constant   4   0    as    k   ■*■    ±   °°. 

l'  r  r 

This  condition  is  represented  in  Figure  13. 


FIGURE  13 
Dispersion  Relation  in  the  Complex  u)  Plane 
for  an  Undefined  Instability 


The  only  conclusion  that  can  be  drawn  in  this  case  is 
that  the  disturbance  is  unstable .   It  is  Sturrock's  conten- 
tion that  this  situation  will  occur  only  when  the  model  is 
inexact  or  incomplete „ 

To  apply  Sturrock's  theory  to  the  WKB  solution  obtained 
for  the  perturbed  wave  function,  let 

r 


k  =  ft  + 


4ft{x2+l)   4ftx 


15 


(x2+l) 


+  ft' 


(16) 
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when  15       . 

ft2(x2+l)2 

and 


/l5  .   -1      ft2 
k  =  tan   x  + 


(4x2+l)  +  i-  In 


2/T5 


M—  +  ft2 


L  (x2+d 


(17) 


when 

ft2(x2+l)2 
15 


<  1   . 


Plots  of  ft.  vs.  ft   and  k.  vs.  k   are  shown  in  Figures  14 
l      r      i      r  ^ 

and  15  for  x  =1.5.   The  plot  of  k. ,  k  was  obtained  directly 
from  equations  (16)  and  (17).   The  plot  of  ft.,  ft  was  obtained 
by  approximating  the  logarithms  and  then  solving  for  ft. 
For  equation  (16)  at  x  =  1 . 5  , 

k  =  ft  +  ii£i  +  i(.166)ln[1.42+ft2]  .  (18) 

The  logarithm  is  approximated  by  factoring  out  ft2  and 
dropping  the  remainder, 

ln[1.42+ft2]  =  21n  ft  +  In  [i-iii  +1]   . 

ft2 

But  from  the  restrictions  on  equation  (18) ,  due  to  aoproxima- 
tions  in  the  last  section,  the  last  term  cannot  be  greater 
than  In  1.63.   Since  for  this  equation  |ft|  >  1.51,  the  last 
term  may  be  dropped  if  ft  is  very  large.   The  logarithm  of  ft 
is  then  approximated  by 

111  "  =  "ft"    * 

Although  this  introduces  a  great  deal  of  error  if  ft  be- 
comes too  large,  it  does  allow  one  to  solve  for  ft  using  only 
a  quadratic  and  the  curve  is  correct  in  general  features. 

The  result  is 

,    ~,    3.61    .  ,  -, -j -> \  (ft- 1) 

k  =  ft  +  — 7=; — ■  +  i(.333)  — q—   / 
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Dispersion  Relation  in  Complex  oo  Plane  for  the 
Perturbation  of  the  Static  Solution 

Figure  14 
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Dispersion  Relation  in  Complex  k  Plane  for  the 
Perturbation  of  the  Static  Solution 

Figure  15 
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which  may  be  solved  for  ft, 

ft  =  ~   -  i(.166)  ±  %\/k2-14.56+|i(2-k)  .     (19) 

The  same  general  ideas  may  be  carried  out  for  equation 
(17)  . 

k  =  2.54  +  .266ft2  +  i(.117)ft2  +  i(.0585)   . 


ft 


=  ±  2.04e~1  11'90°  Vk-2.54-i(.0585)       (20) 


This  equation  is  useful  only  if  |ft|  <  .95.   This  restric- 
tion severely  limits  the  values  allowed  by  k  and  the  equation 
is  not  very  useful  in  obtaining  Figure  14.   The  plot  was  ob- 
tained using  equation  (19).   When  the  ft.,  ft  plot  is  dis- 
placed below  the  ft   axis,  the  motion  is  stable.   This  is  the 
case  for  large  values  of  ft   in  Figure  14.   Nearer  the  origin 
one  can  see  that  the  plot  looks  as  though  it  is  unstable.   In 
this  region  r-  >   1  .   This  corresponds  to  3  >  1  in  section 
2.3.2  where  an  unstable  solution  was  predicted.   The  upper 
curve  in  this  region  has  been  terminated  because  neither  equa- 
tion 19  nor  20  yield  valid  values  of  ft.   However,  the  curve 
seems  to  curve  back  towards  the  origin  from  calculations  with 
equation  (19)  .   This  would  be  consistent  with  the  idea  that 
the  solution  becomes  stable  for  large  values  of  ft. 

The  k.,  k   plot,  Figure  15,  does  not  give  as  much  infor- 
mation.  The  large  gap  between  the  two  line  segments  is  due 
to  the  limitations  on  equations  (19)  and  (20) . 

It  can  be  shown  with  equation  (.16)  that 

,  .   k. 
lim  __i  _  0 

ft*-oo  k 

r 
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However,,  due  to  the  limitations  on  the  ranges  of  equations 
(16)  and  (17) ,  one  cannot  determine  very  much  about  the  aaos 
in  Figures  14  and  15.   Therefore,  one  cannot  say  much  about 
the  type  of  instability  that  exists  for  small  k. 
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CHAPTER  III 

SUMMARY  AND  CONCLUSIONS 

3.1   SUMMARY 

The  nonlinear  wave  equation, 

2 


72$  -  eX^i-K    5LB  o  , 


C   8t2 

derived  from  the  Klein-Gordon  wave  equation  has  been  investi- 
gated.  The  e  =  +1,  this  equation  represents  a  true  analog 
to  the  Klein-Gordon  wave  equation  with  the  mass  replaced  by 
field  self-interaction.   For  e  =  -1,  the  nonlinear  term  may 
be  interpreted  as  the  source  term  in  the  electromagnetic 
potential  equation  with  the  source  given  by  field  self- 
interaction. 

Solutions  to  the  static  wave  equation  are  characterized 
by  a  fundamental  length  k  and  an  "energy  parameter",  h,  as 
integration  constants.   The  term  "energy  parameter"  is  used 
only  in  the  mechanical  analog  sense.   It  is  found  that  for 
e  =  +1,  the  solutions  all  are  unbounded  and  not  everywhere 
real  for  all  values  of  h.   The  case  h  =  0  may  behave  either 
as  a  1/r  potential,  infinite  at  the  fundamental  distance  k, 
and  with  a  nonreal  solution  for  r  <  k ;  or  as  a  potential  well, 
centered  at  the  origin,  and  infinite  at  r  =  K  with  no  real 
solution  for  r  >  k.   The  first  result  may  allow  the  inter- 
pretation of  a  particle  with  a  "hard  core".   The  second  re- 
sult could  be  interpreted  as  the  basis  of  a  particle  whose 
characteristics  are  determined  by  eigenfunctions  within  the 
potential  well. 
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For  e  =  -1,  the  solutions  are  everywhere  real.   In  partic- 
ular „    the  solution  with  h  =  0  is  also  everywhere  finite  and 
has  a  1/r  dependence  for  r  >  k.   This  solution  mav  represent 
an  electromagnetic  model  of  a  charge  with  k  as  the  extension 
of  the  charge. 

Solutions  of  the  full  wave  equation  were  obtained.   The 
solutions  are  distinguished  in  their  character  by  their  veloc- 
ity relative  to  light  and  by  the  sign  of  e.   Those  solutions 
for  which  3  <  1/  £  =  1  or  3  >  1,  e  =  -1,  are  unbounded  and 
not  everywhere  real.   Those  solutions  for  which  6  >  1,  e  =  1, 
or  3  <  1,  e  =  -1  are  periodic. 

The  study  of  the  stability  of  the  solutions  with  z    =  -1 
represented  the  major  objective  of  this  caper.   The  stabilitv 
was  tested  in  three  different  ways.   The  stability  of  the 
periodic  solutions  was  tested  for  one  dimensional  solutions 
using  the  method  of  characteristic  exponents  used  in  non- 
linear mechanics.   These  solutions  were  found  to  be  indif- 
ferent to  perturbations,  that  is  the  perturbation  results  in 
a  new  solution  which  will  be  maintained  without  blowing  up  or 
returning  to  the  original. 

Next,  the  stability  of  a  perturbation  on  the  static  do- 
tential  for  h  =  0  was  investigated.   The  wave  equation  for 
the  perturbation  was  linearized  and  analyzed  for  stability 
assuming  that  cj>Q  is  constant.   This  amounts  to  studying  "local 
stability".   If  B2  <  1,  Liapounoff's  stability  analysis  gives 
no  information  about  the  stability  of  equilibrium.   The  per- 
turbation in  this  case  is  periodic  and  may  be  analyzed  for 
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stability  of  the  focal  point  of  the  phase  diagram  using  the 
method  of  characteristic  exponents.   The  result  found  was 
again  indifferent  to  stability.   If  S2  >  1,  Liapounoff's 
method  indicates  instability. 

Finally,  the  linearized  perturbation  equation  was  inves- 
tigated using  Sturrock's  analysis  of  the  dispersion  relation 
obtained  from  a  WKB  approximation  to  the  solution.   The  dis- 
persion relation  obtained  was  investigated  at  r/<  =  1.5, 
i.e.,  near  the  center  of  the  charge.   Of  the  two  possible 
complex  roots  of  oo  for  real  wave  number,  k,  one  was  found  to 
have  only  negative  imaginary  parts  indicating  a  stable  solu- 
tion.  The  other  solution  has  positive  imaginary  values 
around  go   =0,  indicating  instability  when  6  >  1.   The  char- 
acter of  the  instability  could  not  be  determined  because  of 
restrictions  on  the  equations  used  due  to  the  WKB  and  other 
approximations . 

The  results  obtained  from  Sturrock's  analysis  agrees  with 
the  results  obtained  from  the  first  two  methods,  if  3  >  1. 
For  3  <  1/  Sturrock's  method  predicts  stability,  whereas  the 
method  of  characteristic  exponents  predicts  indifference. 

3.2   CONCLUSIONS 

It  has  been  demonstrated  in  the  paper  that  a  relatively 
simple  nonlinear  field  equation  can  serve  as  a  convenient  test 
ground  for  investigation  of  stability  concerning  solutions  of 
nonlinear  field  equations  which  might  represent  particles  or 
electromagnetic  sources  arising  from  field  self-interaction. 
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The  analysis  in  this  case  has  shown  that  certain  perturba- 
tions are  indifferent  or  stable.   The  result  seems  to  be  sen- 
sitive to  the  change  from  consideration  of  "local  stability" 
to  consideration  of  the  varying  <J)_.   This  suggests  that  a  non- 
linear stability  analysis  might  eventually  change  the  result 
even  further.   To  determine  whether  or  not  this  is  true  will 
require  a  better  approximation  to  the  solution  for  the  oer- 
turbation.   This  may  be  accomplished  by  including  nonlinear 
terms  in  the  linearized  equation  and  by  carry inq  the  WKB 
method  to  high  orders  of  approximation  by  some  kind  of  itera- 
tion.  A  second  alternative  would  be  to  seek  new  avenues  of 
stability  analysis,  eventually  using  computer  techniques. 

It  has  been  shown  that  there  exist  other  classes  of  per- 
turbations which  are  unbounded.   It  remains  to  be  established 
whether  or  not  these  perturbations  are  everywhere  real  and 
what  their  physical  interpretation  would  be,  perhaps  connected 
with  shock-like  phenomenon.   Answers  to  these  questions  aqain 
lie  in  improvement  of  the  WKB  approximation  and  the  analysis 
of  the  dispersion  relation. 
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APPENDIX  A 
The  General  Solution  of  the  Static  Wave  Equation 

In  the  main  body  of  the  paper  ootential  functions  were 
obtained  for  the  case  h  =  0  in  the  equation 


dFj  2  =  2h  +  £^  F   +  %F2 
du 


(4) 


In  general,  the  solutions  of  the  equation  are  elliptic 
functions  whose  form  is  dependent  upon  the  sign  of  h  and  e, 
and  upon  the  relative  size  of  h  and  X. 

In  equation  (4) ,  let   F  =  V 


and 


24h  =  6a2 


where 


6  =  ±1  . 
The  equation  then  becomes 
dvl2 


3(hS)   =  vU4A2v3  +  3v+6a2)  . 


This  equation,  when  it  yields  elliptic  functions,  may  be 
evaluated  using  a  handbook  of  elliptic  functions.  [6]   The 
cubic  on  the  right  must  first  be  factored. 

3    dv  2 


4X2   du 


A+B    (A-B)   /— =■,  r     A+B    (A-B)   /— -,  -, 
ev{v-  (A+B)  ]  [v  +  ^-  -        0    /-3]  [v  +  -*-   +    9    S=3]  ) 


where 


A  "  2X 


6a2X    r 

0  .  +  V-' 


aHXz+e 
o 


B  = 


2X 


6a2X 


-Va^X2  +  £ 


(21) 


(22) 


(23) 
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To  solve  the  equation,  the  solution  shall  be  divided  into  two 
cases . 

Case  1.   e  =  1 
If  6  =  1  (i.e.,  h  >  0)  then  F  can  have  all  values  and 

0  <  v  <  °°  . 
The  differential  equation  is 


2A 


dv 


k 


[v(v-R1) (v-R2) (v-R3) ] 2 


i   r 
In  — 

K 


where  R  has  been  substituted  for  the  roots.   For  this  case, 
A  and  B  are 


h   ["a0X  *Va»X*  +  l]3    . 


From  an  elliptic  functions  handbook  the  solution  to  this  equa- 
tion is  obtained  by  letting 

2    (ot+3)  2-(a+b)2 


kz  - 


cos 


4a3 
2  A  (a- 3)  v-  (a+b)a 


2^ 


* 


Then 


2A(a+3)v-(a+b)a 
a   =  (a2+  1  +  bz) 
3  =  [3(a2  -  1  +  bz)] 
a  =  2AA,  b  =  2AB 


/a3  n   r 
cos  t)  =  en  \/-s—  In  —  , 
3      K 


where  the  modulus  of  the  elliptic  function  is  k. 
Let  the  argument  of  the  en  be  x.   Then 


*0  =  [V 


fl 


\2Ar 

i 


g(a+b)tn2x/2  dn2x/2 
,atn2x/2  dn2x/2-3 


h 
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This  solution  actually  gives  a  series  of  unbounded  waves. 
The  complicated  function  is  real  for  small  values  of 

tn2x/2  dn2x/2 
because  a  +  b  is  negative. 

If  6  =  -1  (i.e.,  h  <  0  the  differential  equation  is 


/5 

2X 
j 


r 


v 

[v(v-R1) (v-R2) (v-R3)]%  =  ln  \ 
Rl 


where  the  roots  are  defined  by  equations  (21) ,  (22) ,  and 

(23) .   Specifically  A  and  B  are 

1 

a   =    2XA  =     [a2X    +  Va"X2+l]     3 


b   =    2XB   =    [a2X    -  VaJX2+l] 

and      R      =   A   +    B 

The  solution  will  be  very  similar  to  that  for  6=1, 
except  that  from  section  1.2.1,  F  has  a  minimum  of 

F  min  =  (A+B) 2 
or  v  min  =  A  +  B. 

This  will  necessitate  redefining  slightly  the  elliptic  func- 
tion and  its  modulus. 

k2  =  _(a_+$)2  -  (a+b)2 


4a6 
a2  =  3(a2-l+b2) 

62  =  a2  +  1  +  b: 


2X (a-B)v  +  (a+b) $ 
cos  9  "  2X(a+g)v  -  (a+b) 3 
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fa  6  r 

cos    0   =   cn\/-4  In  —  =  en   x 
y     3  k 


^0  =  2X7 


B(a+b) 


3-atn2x/2  dn2x/2/_ 


This  function  also  represents  an  unbounded  wave.   The  wave 
minimum  occurs  periodically  and  approaches  zero  as  r  increases 

Case  2.   e  =  -1 

If  6  =  1,  the  differential  equation  will  be 


7  ia^)2  =  -v[(v-R1)(v-R2)(v-R3)] 


4A 

where  R, ,  R-»  R~  are  defined  by  equations  (21) ,  (22) ,  and 

(23) .   A  and  B  will  be  , 

a  =  2AA  =  [a2  A  +  ^fa^XT^l]     3   , 

1 
b  =  2AB  =  [a2A  -Va^A"2™]  3   . 

In  this  instance  there  are  three  different  solutions.   Assume 
that 

ajA2  >  1   . 


Then 


/3 

2A 


v 


dv ,    ,   r 

[v(R,-v)  (v-R~)  (v-RJ]  *,  ~    ln  k    , 


where  the  roots  will  be  specified  as 


R,  =  A  +  B   and 


R„  and  R_.  are  the  complex  roots. 
From  section  1.2.2,  V  has  the  range 

0  <•  v  <  A  +  B  . 
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Let 


k2  = 


(a+b)2-(a-B)2 

£aB 


cos  6  = 


(a+b)  B-2Av(B+a) 
(a+b) B+2Av(a-B) ' 


Then 


a2  =  3(a2+l+b2),  B2  =  a2  -  1  +  b: 


cos  6  =  cnx/^4  In  -  =  en  x  , 

y  3       k 


where  the  elliptic  function  has  the  modulus  k.   From  the  last 
equation  one  may  solve  for  v  and  cj>0, 


*0  " 


1   (B (a+b) tn2x/2dn2x/2' 
2Xr  U+Btn2x/2dn2x/2 


h 


<j>0  oscillates  in  this  case  with  decreasing  amplitude  as  r  in- 
creases . 

Assume  a)U2  =  1-  This  is  a  special  case  in  which  the 
solution  may  be  expressed  in  terms  of  elementary  functions. 
One  has 


/3 
2A 


r 


0  [v(f-v)(v+^)2]' 


dv         ,   r 
— i r  =  In  —  . 

lx  2-.%  < 


Make  the  substitution 

1 

v  =  y  -  it  • 

Then  from  a  table  of  integrals: 


/3 
2X 


c- 


*L 


%X  Yl-4X2+  X 


2v_,-2  h 


-v 


sin 


4Ay-3 
2Xy 


Y 
J  %X 


.  -1 


Upon  transposing  the  sin    and  some  manipulation, 


* 


0 


1    1-cos  In  r/K 
2Xr   2+cos  In  r/K 
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Assume   an^2  >    !•   A  better  representation  of  the  roots 
of  the  cubic  in  this  case  would  be 


R. 


y-  cos  cj>/3 


R2   =    — ~y    (cos    <J)/3+/3sin   cf>/3)       , 
R3   =      2T    ("cos    ^/3+/3sin    <J>/3)     , 
cos    (j)    =   a^A    . 

The    range   of   v   as    calculated   in    section    1.2.2    is 

1 

0    <_  v  _<   t-  cos    (J)/3    . 


The   differential   equation   is 
v 


/3    f  dy_ 

2XJ  1/X    cos    */3    [v(R1-v)(v-R2)(v-R3)]- 


i       r 
r-     =    In  — 


Define 


k2    = 


n  ■  2(f) 

2    sin  — r 


/?  2<j>  2$ 

/3    cos    — ■=■  +    sin  — ^ 


sin 


/[cos    (j)/3+/3sin   4>/3]  [cos    cfr/3-Av] 
V      cos    <J>/3  [2Av+cos    (f>/3+/3sin    4>/3] 


Then 


or 


sin    9    =    sn   % 


/3cos 

20/3+sin 

2^/3] 

h 

In 

K 

i 

/3 

J 

r 

sin 

0    =    sn 

X 

• 

From   the    last   equation   v   and    <ju   may   be   obtained 


_1 

Ar 


en2  Acqs    cj)/3(cos    cj)/3+/3'sin    4>/3) 
.cos    4>/3(l+2sn2A)+/3~sin    <t>/3 


% 


Let  6  =  -1.   Then  the  differential  equation  is 


3 

4A2 

rdvl 
LduJ 

=  -V 


3v    "0 

4T2    xr2 


IT 
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A  and  B  in  the  roots  as  defined  by  equations  (21) ,  (22) ,  and 


(23)  are 


a  =  2XA  =  [-a2X  +  "Va^X2-!] 


b  =  2XB  =  [-a2X  -  ^JakQXz-l) 


1 
3 

1 
I 


It  was  shown  in  section  1 „ 2 „ 2  that  when  h  <  0,  F  cannot 
be  real  unless,   aQ^2  .1  ^°   Therefore  this  will  be  the  only 
case  considered.   The  roots  can  be  better  expressed  as 

R,  =   t-  cos  <t>/3  , 

R-,   =    -Tj-    (-cos    <J)/3-/3sin    c{>/3)       , 


2X 


R. 


Yr    (-cos    4>/3+/3sin   <J>/3) 


and   v  has    the    range 


cos    <$>   =   ~a~.X    , 


R3    <    v    <    L 


Therefore 


/3 
2X 


v 


R. 


dv 
[v(R,-v)  (v-RJ  (v-R-J  ] 


%   =   In   - 


Define : 


2         /3cos    20/3+sin    2<j>/3 
K  2sin    2<t>/3 


sin    6    =- 


4/3sin    <f>/3(cos    <j>/3-Xv) 


V   (3cos    <j»/3-/3sin   <(>/3)  (2Xv+cos    <{>/3+/3sin    <(>/3) 


Then 


sin    8    =   sn 


sin    24>/3l  *    ,       k 

— !s  ln  r 


=   sn  x 
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^0  =  mr 


where  the  sn  x  has  modulus  k.   By  squaring  both  sides  and 
solving  for  v  one  may  obtain 

1   f  sin  2<j)/3(2-sn2x)-/3sn2xcos  2<1>/3U  * 
sin  4)/3  (2-sin2x) +/3sn2xcos  (J)/3  /I 

This  is  an  oscillating  function.   In  the  limit  as  cos  9  =  -1, 
the  term  in  the  brackets  becomes  identially  1.   Then 

1 


*0  = 


i  % 


2Xr 


and  F  is  a  constant  as  predicted  in  section  1.2.2. 
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APPENDIX  B 

Exact  Solution  of  /Kdx 

The  exact  solution  of  this  integral  involves  all  three 
types  of  elliptic  functions.   One  may  obtain  the  solution  in 
the  following  manner. 


k 


Kdx  = 


[15+ft2  (x"  +  ir  ]  2dx 
x2  +  l 


Let   z  =  x2  +  1.   Then  the  integral  becomes 


[15+ft2z2]  %dz 


2z(z-l) 


k 


This  may  be  reduced  to 


15dz 


+ 


fiz  dz 


2fiz(z-l)%  l±l   +  z2)%     2(z-l)%/15 


+  z 


\T 


n* 


Define : 


3  =  (15+fi2)%  , 


k2  =  3  "n 


cos  0  = 


23 

3  +  fl-gz    3  -  ftx2 
3  "  "+Uz  "  3  +  «x2 


v  =  F(9,k) 
where  F(9,k)  is  the  elliptic  integral  of  the  first  kind, 
The  integrand  may  be  replaced  by 


k 


(1+cn  v)dv      dv  ,Q,  ^     [3+ft+(ft-3)cn  v] 


11    (SL\ 
2     V   3+ft+(ft-B)cn  v 


+  -^  (F) 


1+cn  v 


where  en  v  is  an  elliptic  function  of  v. 
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Integrating  over  v  one  obtains 


v 


ru%  15 


3  '      B-fi 


10    n\%i-./    \  /do\%   sn   vdn   v 

-    (3-ft)    E(v)    +    (3ft)      1+cn     — 


15 


4(ft3) 


3+ft 
%     13-ft 


.(9,    "^    ,    k) 


3       .       -1  f  15      2 

+  ~  tan       !    43      M 


sd   v 


where   E(v)    is    the    elliptic   integral   of   the   second  kind, 

-3k' 


7T(9, 


n 


,k)  is  the  elliptic  integral  of  the  third  kind,  and 


the  functions  sn  v,  dn  v,  sd  v,  are  elliptic  functions  of  v. 

The  following  equalities  apply: 

sn  v  =  sin  0  , 

en  v  =  cos  G  , 

dn  v  =  (l-k2sin29)^, 
sin  9 


sd  v  = 


k 


(l-k2sin29)^ 

Using  these  equalities  the  elliptic  functions  may  be  replaced 
by  functions  of  x.   The  results  yield: 


v  «f)  %  F^T"  ^ft)W)  - 


15 


4(ns) 2 


3+ft 
3-ft 


ir(e,-S£l,k) 


+_§.  tan"1 
/2 


15  x 


k 


.223[32+2ft2x2+ft2x't] 


+ 


ft 


3+ftx' 


(32  +  2Q2x2+ft2x'+) 


k 
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